rm(list = ls())
library(spatialreg)
library(cshapes)
library(tidyverse)
library(sf)
library(stargazer)
library(lubridate)
library(tscsdep)
library(texreg)
library(xtable)
library(texreg)
library(plm)
library(spdep)

load("data/data_imputed.RData")
source("code/st_effects_new.R")


data  <- data_imputed
data <- data %>% 
       dplyr::rename(country = country_name) 

## lag all IV for one year
data <- data %>% 
       dplyr::arrange(cowcode, year) %>% 
       dplyr::group_by(cowcode) %>% 
       dplyr::mutate_at(vars(c("logpop", "loggdppc", "tradependonCHN",
                               "v2x_polyarchy", "polity2", "capratio",
                               "dip_hist", "USpresident_visit","solschange",
                               "statehead.total", "statefrom.total",
                               "stateto.total")), ~dplyr::lag(.,n = 1L))  

data <- data %>% 
       dplyr::arrange(cowcode, year) %>% 
       dplyr::group_by(cowcode) %>% 
       dplyr::mutate(laglogidealdistance = dplyr::lag(logidealdistance, n = 1L),
                     lagagree = dplyr::lag(agree, n=1L))
       
data <- data %>% 
       dplyr::filter(!is.na(statefrom.total)) %>% 
              ungroup()

########################## Table 1 ##########################

# DV: Distance in Ideal Points
# Model 1; OLS
reg1 <-lm(formula = logidealdistance ~ laglogidealdistance + statefrom.total +
                logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                dip_hist + solschange + factor(cowcode) + factor(year), 
         data = data)

## Model 2: SAR
wm1 <- make_ntspmat(lmobj = reg1, ci = "country", yi = "year", k = 5)
sar1 <- ntspreg(reg1, wm = wm1)
sac1 <- ntspsac(lmobj= reg1, wm = wm1) 
err1 <- ntsperr(lmobj = reg1, wm = wm1)
## Model 3: SLX
weights1 <- mat2listw(wm1[[2]], row.names = NULL, style = "W")
nblist1 <- weights1$neighbours
listw1 <- nb2listw(nblist1)
slx <- lmSLX(logidealdistance ~ laglogidealdistance + statefrom.total +
                    logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                    dip_hist + solschange + factor(cowcode) + factor(year), data = data,
             listw = listw1, Durbin = ~statefrom.total)


### DV: Agreement with China
# Model 4; OLS
reg2 <-lm(formula = agree ~ lagagree + statefrom.total +
                 logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                 dip_hist + solschange + factor(cowcode) + factor(year), 
          data = data)

# Model 5; SAR
wm2 <- make_ntspmat(lmobj = reg2, ci = "country", yi = "year", k = 5)
sar2 <- ntspreg(lmobj = reg2, wm = wm2)
sac2 <- ntspsac(lmobj= reg2, wm = wm2) 
err2 <- ntsperr(lmobj = reg2, wm = wm2)
# Model 6; Slx
weights2 <- mat2listw(wm2[[2]], row.names = NULL, style = "W")
nblist2 <- weights2$neighbours
listw2 <- nb2listw(nblist2)

slx2 <- lmSLX(agree ~ lagagree + statefrom.total +
                     logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                     dip_hist + solschange + factor(cowcode) + factor(year), data = data,
              listw = listw2, Durbin = ~statefrom.total)

############## Output to Table 1 ############## ##############

texreg(list(reg1, sar1, slx, reg2, sar2, slx2), file = "tables/mod1.tex",
       stars = c(0.01, 0.05, 0.1),
       caption = "The Impact of Chinese Leadership Visits on Foreign Policy Convergence",
       caption.above = TRUE,
       label = "tab:tscsm1",
       scalebox = 0.8,
       custom.header = list("DV: Distance in Ideal Points" = 1:3,
                            "DV: Agreement with China" = 4:6),
       custom.model.names = c("M1:OLS", "M2:SAR", "M3:SLX", "M4:OLS", "M5:SAR","M6:SLX"),
       custom.coef.map = list("laglogidealdistance" = "Temporal Lagged DV",
                              "lagagree" = "Temporal Lagged DV",
                                 "statefrom.total" = "Chinese State Visits",
                              "$\\rho$" = "Spatial Lag Parameter: $\\rho$",
                             # "$\\lambda$" = "Spatial Error Parameter: $\\lambda$",
                              "lag.statefrom.total" = "$\\boldsymbol{W_x}$: Chinese State Visits",
                                 "logpop" = "Ln(Population)",
                                 "loggdppc" = "Ln(GDP per Capita)",
                                 "tradependonCHN" = "Trade Dependence",
                                 "v2x_polyarchy" = "Electoral Democracy Index",
                                 "capratio" = "Capability Ratio",
                                 "dip_hist" = "Years of Diplomatic Relations",
                                 "solschange" = "Change in Source of Leader Support"),
          custom.gof.rows = list("Fixed Country Effects" = c("YES", "YES", "YES", "YES", "YES", "YES"),
                                 "Fixed Year Effects" = c("YES", "YES", "YES", "YES", "YES", "YES")),
          digits = 4)

############## output to Table 2 ############## ##############
# Table 2. The Short-run and Long-run Effects of Chinese Leadership Visits

#st_effects calculates both the short and long-run total, direct and indirect effects of shocks to sample units' covariates using the last cross-section of the data.
eff <- st_effects_new(spobj = sar1, wm = wm1, tlag = "laglogidealdistance", covariate = "statefrom.total")
print(xtable(eff, caption = "The Short-run and Long-run Effects of Chinese Leadership Visits",
             label = "tab:eff"),include.rownames=FALSE,
      file = "tables/sareffect.tex")
########################################################



############## Figure 5. The Robust Effects of Chinese Leadership Visits using K =2, · · ·, 20

ks = seq(2, 20, by = 1)
wms <- list()
for (i in seq_along(ks)){
       wms[[i]] <- make_ntspmat(lmobj = reg2, ci = "country", yi = "year", k = ks[i])     
}


library(purrr)
wmsk <- list()
for (i in seq_along(ks)){
       wmsk[[i]] <- wms[[i]][[2]]       
}

## modeling
slxs <- list()

for (i in seq_along(ks)){
 
weights_k <- mat2listw(wmsk[[i]], row.names = NULL, style = "W")
nblist_k <- weights_k$neighbours
listw_k <- nb2listw(nblist_k)      
       
slxs[[i]] <- lmSLX(logidealdistance ~ laglogidealdistance + statefrom.total +
                    logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                    dip_hist + solschange + factor(cowcode) + factor(year), data = data,
             listw = listw_k, Durbin = ~statefrom.total)
}


kb <- slx_plot_knb(modResults = slxs, var = "lag.statefrom.total") + 
       xlab("K") + ylab("Wx: Chinese State Visits") 
ggsave("kn.pdf", height =6, width = 8)

kb2 <- slx_plot_knb(modResults = slxs, var = "statefrom.total") + 
       xlab("K") + ylab("Chinese State Visits") 
ggsave("kn2.pdf", height =6, width = 8)



### Appendix: Table A1, A2, A3, Figure A1

### Table A1. Descriptive Statistics of Variables
## get a summary table
sum_df <- data %>%dplyr::select(logidealdistance,
                             agree,
                             statefrom.total,
                             logpop, loggdppc,
                             tradependonCHN,v2x_polyarchy,
                             capratio,
                             dip_hist, solschange, 
                             polity2,USpresident_visit)
sum_df <- as.data.frame(sum_df)
library(stargazer)
stargazer(sum_df, type = "latex", title = "Descriptive Statistics of Variables",
          digits = 4, header = FALSE, label = "tab:sum"
          )


#Table A2. Robustness Check: The Impact of Chinese Leadership Visits on Foreign Policy Convergence

texreg(list(err1,sac1,err2, sac2), file = "tables/mod2.tex",
       stars = c(0.01, 0.05, 0.1),
       caption = "Robustness Check: The Impact of Chinese Leadership Visits on Foreign Policy Convergence",
       caption.above = TRUE,
       label = "tab:tscsm2",
       scalebox = 0.9,
       custom.header = list("DV: Distance in Ideal Points" = 1:2,
                            "DV: Agreement with China" = 3:4),
       custom.model.names = c("Model1:SEM", "Model2:SAC", "Model3:SEM", "Model4:SAC"),
          custom.coef.map = list("laglogidealdistance" = "Temporal Lagged DV",
                                 "lagagree" = "Temporal Lagged DV",
                                 "statefrom.total" = "Chinese State Visits",
                                 "$\\rho$" = "Spatial Lag Parameter: $\\rho$",
                                 "$\\lambda$" = "Spatial Error Parameter: $\\lambda$",
                                 "lag.statefrom.total" = "$\\boldsymbol{W_x}$: Chinese State Visits",
                                 "logpop" = "Ln(Population)",
                                 "loggdppc" = "Ln(GDP per Capita)",
                                 "tradependonCHN" = "Trade Dependence",
                                 "v2x_polyarchy" = "Electoral Democracy Index",
                                 "capratio" = "Capability Ratio",
                                 "dip_hist" = "Years of Diplomatic Relations",
                                 "solschange" = "Change in Source of Leader Support"),
          custom.gof.rows = list("Fixed Country Effects" = c("YES", "YES", "YES", "YES"),
                                 "Fixed Year Effects" = c("YES", "YES", "YES", "YES")),
          digits = 4)

## Table A3. Robustness Check: The Impact of Chinese Leadership Visits on Foreign Policy Convergence

ols_rob1 <-lm(formula = logidealdistance ~ laglogidealdistance + statefrom.total +USpresident_visit+
                 logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                 dip_hist + solschange + factor(cowcode) + factor(year), 
          data = data)

wm1_rob <- make_ntspmat(lmobj = ols_rob1, ci = "country", yi = "year", k = 5)

sar1_rob1 <- ntspreg(reg1, wm = wm1_rob)

weights_rob1 <- mat2listw(wm1_rob[[2]], row.names = NULL, style = "W")
nblist_rob1 <- weights_rob1$neighbours
listw_rob1 <- nb2listw(nblist_rob1)

slx_rob1 <- lmSLX(logidealdistance ~ laglogidealdistance + statefrom.total +USpresident_visit+
                         logpop + loggdppc + tradependonCHN + v2x_polyarchy + capratio + 
                         dip_hist + solschange + factor(cowcode) + factor(year), data = data,
              listw = listw_rob1, Durbin = ~statefrom.total) 

###### 
ols_rob2 <-lm(formula = logidealdistance ~ laglogidealdistance + statefrom.total +
                     logpop + loggdppc + tradependonCHN + polity2 + capratio + 
                     dip_hist + solschange + factor(cowcode) + factor(year), 
              data = data)

wm2_rob <- make_ntspmat(lmobj = ols_rob2, ci = "country", yi = "year", k = 5)

sar1_rob2 <- ntspreg(ols_rob2, wm = wm2_rob)

weights_rob2 <- mat2listw(wm2_rob[[2]], row.names = NULL, style = "W")
nblist_rob2 <- weights_rob2$neighbours
listw_rob2 <- nb2listw(nblist_rob2)
slx_rob2 <- lmSLX(logidealdistance ~ laglogidealdistance + statefrom.total +
                         logpop + loggdppc + tradependonCHN + polity2 + capratio + 
                         dip_hist + solschange + factor(cowcode) + factor(year), data = data,
                  listw = listw_rob2, Durbin = ~statefrom.total) 

texreg(list(ols_rob1,sar1_rob1, slx_rob1, ols_rob2,sar1_rob2, slx_rob2),
       file = "tables/mod4.tex",
       stars = c(0.01, 0.05, 0.1),
       caption = "Robustness Check: The Impact of Chinese Leadership Visits on Foreign Policy Convergence",
       caption.above = TRUE,
       label = "tab:tscsm4",
       scalebox = 0.8,
       custom.header = list("DV: Distance in Ideal Points" = 1:6),
       custom.model.names = c("Model1:OLS","Model1:SAR", "Model3:SLX", "Model4:OLS", "Model5:SAR",  "Model6:SLX"),
       custom.coef.map = list("laglogidealdistance" = "Temporal Lagged DV",
                              "statefrom.total" = "Chinese State Visits",
                              "$\\rho$" = "Spatial Lag Parameter: $\\rho$",
                              "lag.statefrom.total" = "$\\boldsymbol{W_x}$: Chinese State Visits",
                              "USpresident_visit" = "US Presidential Visits",
                              "logpop" = "Ln(Population)",
                              "loggdppc" = "Ln(GDP per Capita)",
                              "tradependonCHN" = "Trade Dependence",
                              "v2x_polyarchy" = "Electoral Democracy Index",
                              "polity2" = "Polity IV Scores",
                              "capratio" = "Capability Ratio",
                              "dip_hist" = "Years of Diplomatic Relations",
                              "solschange" = "Change in Source of Leader Support"),
       custom.gof.rows = list("Fixed Country Effects" = c("YES", "YES", "YES", "YES","YES", "YES"),
                              "Fixed Year Effects" = c("YES", "YES", "YES", "YES","YES", "YES")),
       digits = 4)



## Figure A1 shows the correlation coefficients for the covariates.
library(GGally)
ggpairs(sum_df)
library(ggcorrplot)
p <- sum_df %>%  cor(use = "pairwise") %>%
       round(2) %>% 
       ggcorrplot(., type = "lower", lab = T, show.legend = T,
                  legend.title = "correlation\n coefficient",
                  sig.level = 0.05) 
   p +  scale_x_discrete(labels = c("agreement score", "state visits", "population",
                                    "GDPpc","trade dependence", "electoral demo index", 
                                    "capability ratio", "diplomatic relations",
                                    "SOLS change", "Polity IV","US visit"))+ 
          scale_y_discrete(labels = c("ideal point", "agreement score", "state visits", "population",
                                      "GDPpc","trade dependence", "electoral demo index", 
                                      "capability ratio", "diplomatic relations",
                                      "SOLS change", "Polity IV"))
ggsave("cor.pdf", width = 8, height = 8)          



